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Abstract. The gas dynamics in the Galactic disc is modeled by 
releasing an initially axisymmetric SPH component in a com- 
pletely self-consistent and symmetry-free 3D A^-body simula- 
tion of the Milky Way in which the stellar components display 
a COBE-like bar. The density centre of the stellar bar wanders 
around the centre of mass and the resulting gas flow is asym- 
metric and non-stationary, reproducing the HI and CO i — V 
diagrams only at specific times and thus suggesting a transient 
nature of the observed inner gas kinematics. 

The best matching models allow a new and coherent in- 
terpretation of the main features standing out of the i — V 
data within the bar region. In particular, the i — V traces of 
the prominent offset dustlanes leading the bar major axis in 
early-type barred spirals can be unambiguously identified, and 
the 3-kpc arm and its non-symmetric galactocentric opposite 
counterarm are the inner prolongations of disc spiral arms pass- 
ing round the bar and joining the dustlanes at very different 
galactocentric distances. Bania's clump 1 and 2, and another 
velocity-elongated feature near i = 5.5°, are interpreted as 
gas lumps crossing the dustlane shocks. The terminal velocity 
peaks near ^ = ±2.5° are produced by gas along the dustlanes 
and not by the trace of the cusped xi orbit, which passes far- 
ther away from the Galactic centre. According to these models 
and to related geometrical constraints, the Galactic bar must 
have an inclination angle of 25° ± 4°, a corotation radius of 
4.0 — 4.5 kpc and a face-on axis ratio b/a ^ 0.6. 

Key words: Galaxy: structure - Galaxy/ISM: kinematics and 
dynamics - Galaxy: centre - Hydrodynamics - Methods: nu- 
merical 



1. Introduction 

The gas kinematics near the Galactic plane, mainly observed 
in HI (Burton 1985; Kerr et al. 1986; Stark et al. 1992 and 
Hartmann & Burton 1997), ^^CO (Dame et al. 1987; Oka et al. 
1998b), i^CO and CS (Bally et al. 1987), is probably among 
the best tracer of the dynamical mass in the Galaxy. Under the 
assumption of circular motion in an axisymmetric potential, the 
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tangent point method allows in principle to derive the rotation 
curve inside the solar circle and to recover the spatial location 
of the spiral arms (see the compilation of Vallee 1995). 

However, this approach is limited for at least two reasons 
related to the data themselves. First, it is long known (Rougoor 
& Oort 1960) that the longitude- velocity (£ — V) distribution 
of the Galactic gas reveals substantial diffuse and feature-like 
emission in the regions > 0, V < 0) and {£ < 0,V > 0) 
close to the Galactic centre which are forbidden to pure circular 
motion, the most popular example being the 3-kpc arm. Sec- 
ond, bumps of about 7 kms~^ on the terminal velocity curves 
near the tangent points of the spiral arms reflect the gravita- 
tional perturbation of these arms on the gas flow and can prop- 
agate into distance errors of ~ 1 kpc along the line of sight if 
modeled by circular motion (Burton 1992; Combes 1991). 

Traditional interpretations of the "forbidden" velocities 
near the centre involve explosion induced expansion (van der 
Kruit 1971; van der Kruit et al. 1972; Cohen 1975; Oort 1977) 
or gas moving on elliptical orbits in either a spiral (Shane 1972; 
Simonson & Mader 1973) or a barred potential (Peters 1975; 
Cohen & Few 1976; Liszt & Burton 1980; Gerhard & Vi- 
etri 1986). In the last decade, the latter interpretation received 
strong support by direct evidence of a large scale stellar bar 
in the Milky Way from near-IR surface photometry (Blitz & 
Spergel 1991; Dwek et al. 1995; Binney et al. 1997), discrete 
source counts (Nakada et al. 1991 and Izumiura et al. 1994 
for SiO masers, Weinberg 1992 and Nikolaev & Weinberg 
1997 for AGB stars; Whitelock & Catchpole 1992 for Miras 
variables; Blitz 1993 for globular clusters; Sevenster 1996 for 
OH/IR stars; Stanek et al. 1997 for red clump stars) and large 
microlensing optical depths towards the bulge (Paczynski et 
al. 1994; Evans 1994; Zhao et al. 1996; Zhao & Mao 1996; 
Gyuk 1999; see Gerhard 1996 and Kuijken 1996 for reviews). 

Many hydrodynamical simulations (Roberts et al. 1979; 
van Albada 1985; Mulder & Liem 1986; Athanassoula 1992; 
Englmaier & Gerhard 1997) have shown that the gas flow in 
a rapidly rotating barred potential is driven by strong shocks 
leading the bar major axis and followed downstream by en- 
hanced gas densities. In external early-type (i.e. Sbc and ear- 
lier) barred spirals, these shocks are detected as offset dustlanes 
with very large velocity changes in the associated gas velocity 
field (e.g. Reynaud & Downes 1998; Laine et al. 1999). 
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Binney et al. (1991) have interpreted the gas kinematics 
near the Galactic centre approximating the gas streamlines by 
a sequence of closed xi and X2 orbits (Contopoulos & Mertz- 
anides 1977), but such a model does not properly cares about 
the shocks. The hydro simulations designed for the Milky Way 
(Mulder & Liem 1986; Wada et al. 1994; Englmaier & Gerhard 
1999; Weiner & Sellwood 1999) always assume rigid barred 
potentials and bisymmetry with respect to the Galactic centre, 
but some details in the HI and CO observations clearly betray 
important asymmetries: (i) about 3/4 of the molecular gas in 
the nuclear ring/disc lies at positive longitude, (ii) the 3-kpc 
arm and its far-side counterarm have very different absolute 
velocities at ^ = 0, and (iii) the positive and negative veloc- 
ity peaks of bright emission, at {i,V/kms~^) ^ (3°, 270) 
and (—2°, —220) respectively, differ by roughly 50 kms~^ in 
absolute velocity, and the emission of the nuclear ring/disc 
is also shifted towards receding velocities. Moreover, as an 
SBbc galaxy (Sackett 1997), the Milky Way is also expected 
to present prominent dustlanes, but their gaseous traces have 
never been convincingly isolated in the HI and CO data. 

This paper follows a first paper (Fux 1997, hereafter pa- 
per I) where several A/'-body barred models of the Galaxy have 
been built from bar unstable axisymmetric initial conditions, 
constraining the position of the observer relative to the bar with 
the CORE K-band data corrected for extinction by dust. The 
best matching models suggested an angle of 28° ± 7° for the 
bar inclination and a corotation radius of 4.3 ± 0.5 kpc. Here 
we take advantage of these simulations to seek the most conve- 
nient initial conditions for much larger simulations, including 
a gas component treated by the smooth particle hydrodynamics 
(SPH) technique and with no imposed symmetries. The result- 
ing self-consistent gas flows are used to interpret the observed 
gas kinematics near the Galactic plane, and more especially 
within the bar region. Contrary to Weiner & Sellwood (1999), 
our approach is based on the observed bright £ — V features, 
and not on the mostly low density gas close to the terminal ve- 
locities, for which Eulerian codes are more appropriate. 

The structure of the paper is organised as follows: in Sect. 2 
we describe the main features appearing in the HI and CO i — V 
diagrams within the bar region. In Sect. 3 we expose the nu- 
merical A^-body and SPH code used to perform our simula- 
tions. In Sect. 4 and 5, we present the initial conditions and 
the time evolution of these simulations. In Sect. 6 we select 
some optimum models from the simulations and interpret the 
observed gas kinematical features. In Sect. 7 we set some new 
constraints on the bar parameters based on this interpretation. 
Finally, Sect. 8 summarises our results. Throughout the paper, 
we will adopt a solar galactocentric distance Ro = S kpc. 

2. Main observed features 

In the region \£\ < 30° several features stand out from the 
observed longitude-velocity distribution of both atomic and 
molecular gas (see Figs. |l] and||). The nomenclature of these 
features principally refers to Rougoor (1964) and we omit here 
the historical but confusing "expanding" qualification, since lo- 



cal radial motions do not necessarily imply a net flux when 
averaged over azimuth. Other more exhaustive inventories of 
observed i — V features can be found in van der Kruit (1970) 
and Cohen (1975) for the HI, and in Bania (1977), Bally et 
al. (1988) and the review of Combes (1991) for the CO. 

- The 3-kpc arm. This arm is the largest apparent feature, 
covering more than 35° in Galactic longitude, and with a 
"forbidden" radial velocity of —53 km s~^ at ^ = 0. It was 
discovered by van Woerden et al. (1957) and its name is 
related to the location of its tangent point near I = —22°, 
corresponding to R = 3 kpc (at this time the galactocen- 
tric distance of the Sun was assumed to be 8.2 kpc). The 
same authors also noticed that this arm must lie in front 
of the Galactic centre because of absorbing the continuum 
spectrum of radio sources close to the centre (see the HI 
i — V diagram at 6 = in Fig. which is consistent with 
its rather large latitudinal width. The 3-kpc arm cannot be 
represented by a uniformly rotating and expanding circular 
arc over its whole longitude range (Burke & Tuve 1964). It 
was successfully modeled by Mulder & Liem (1986) as a 
stationary density wave in a rotating barred potential. 

- The 135-kms~^ arm. This feature shows no absorption 
lines against Galactic centre radio sources and is sus- 
pected to be the far- side counterpart of the 3-kpc arm (e.g. 
Oort 1977). However it crosses the zero longitude axis at 
135 km s~^ (hence his name), i.e. more than twice the abso- 
lute velocity of the 3-kpc arm, which is among the strongest 
evidence for an asymmetric spiral structure in the central 
few kpc (velocity asymmetries at ^ = cannot result from 
perspective effects). The 135-kms~^ arm extends down 
to ^ ^ —5°, where its radial velocity still reaches about 
100 kms~^, and lies slightly above the Galactic plane, at 
h ^ 0.5°. It also involves less HI mass than the 3-kpc arm 
and is more lumpy. 

- The connecting arm. This is a very rarely discussed fea- 
ture in the recent literature, despite its substantial bright- 
ness. It was probably mentioned the first time by Rougoor 
& Oort (1960). The origin of its name comes from the fact 
that it seems to link the nuclear ring/disc to the outer disc. 
It becomes easily identifiable in the HI and ^^CO data near 
i = 10° and V = 110 km s~^, where it is located far below 
the plane at 6 ^ —1°, and passes through the peak of the 
positive terminal velocity curve. In Rougoor's (1964) de- 
scription, this feature is doubtfully combined with another 
distinct feature extending at roughly constant velocity out 
to ^ ^ 18°. The connecting arm was soon recognized as 
a very inclined spiral arm, i.e. with a high pitch angle, al- 
though its situation in front or behind the Galactic centre 
remained unclear (Rougoor 1964; Cohen & Davies 1976). 
Kerr (1967) adopted the latter alternative and interpreted 
this arm as part of a central gaseous bar. 

- The nuclear ring/disc, or central molecular zone. Con- 
cerns the offcentred dense molecular gas located within 
— 1° ^ i ^ 1.5°, i.e. somewhat inside the positive and neg- 
ative peaks of the terminal velocity curves (see for instance 
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Fig. 1. Grey scale longitude- velocity diagram of ^^CO J = 
1^0 emission within —30° < i < 35° and over a 4° latitude 
strip centred on the Galactic plane (Dame et al. 1999), with the 
dominant features indicated by dashed lines and Bania's (1977) 
molecular clumps. 

Fig. 2. (right column) Longitude- velocity diagrams of HI 
21 cm (left) and ^^CO J = 1 ^ (right) emission as a func- 
tion of latitude in the range \b\ < 1.5°. The HI data are from 
Hartmann & Burton (1997), Burton & Liszt (1978) and Kerr et 
al. (1986), and the CO data from Dame et al. (1987). 

the high resolution CO £ — V maps in Oort 1977; Liszt 
& Burton 1978; Bally et al. 1988; Combes 1991; Morris 
& Serabyn 1996 and Oka et al. 1998b). The parallelogram 
bounding this structure is also called "expanding molecu- 
lar ring" owing to its original interpretation (Scoville 1972; 
Kaifu et al. 1972), and more recently "180-pc ring" accord- 
ing to its size. The molecular complex enclosed by this ring 
forms a kind of disc rotating with a maximum velocity of 
order 100 km s~^ and is speculated to form two mini spirals 
(Sofue 1995). Almost the entire CS emission comes from 
this disc. Binney et al. (1991) have interpreted the paral- 
lelogram and the disc respectively as gas on the cusped xi 
orbit and on X2 orbits in the bar potential (see Sect. |63| ). 
An extensive review on this pattern is given by Morris & 
Serabyn (1996). 
- The molecular ring. Meant to represent a ring-like gas con- 
centration with a radius of roughly 4 kpc, the spatial dis- 
tribution of this structure is in fact poorly known and could 
as well involve imbricated spiral arms. Close to the positive 
terminal velocity curve, it forms two branches with tangent 
points at ^ ^ 25° and 30°, corresponding to ^ 3.4 kpc 
and 4 kpc. The molecular ring probably compares to the 
inner (pseudo-) rings seen in external spiral galaxies (Buta 
1996). 

All features listed here appear both in the CO and HI data with 
a rigorous coincidence in position, although the central molec- 
ular zone is much weaker in HI. It is not yet clear whether 
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these features are transient or represent a permanent gas flow 
on closed orbits. 

Bania (1977) has also isolated two particular "clumps" in 
the ^^CO data (see Fig. |l]), corresponding to massive molecular 
cloud complexes. 

- Clump 1 is located at the negative longitude end of the 
135-kms~^ arm, dXl ^ —5° and h ^ 0.4°, and represents 
the most extreme case of non-circular motion in the Milky 
Way's dense gas kinematics. It is located at slightly lower 
latitude than the nearby emission from the 135-km s~^ arm. 
The clump 1 complex has been resolved into three distinct 
sub-complexes each of about 2.5 x 10^ Mq H2 mass. More 
details on its physical properties are presented in Bania et 
al. (1986). 

- Clump 2 is confined near ^ = 3° and h = 0.2° and spans 
an extraordinary large velocity range of over 150 km s~^ . It 
has been resolved into 16 emitting cores of ~ 5 x 10^ Mq 
by Stark & Bania (1986), who also interpret this clump as a 
dustlane or an inner spiral arm in a barred potential formed 
by molecular clouds distributed along the line of sight. 

The CO is usually considered as a tracer of H2, although it 
is now well established that the density ratio between these two 
molecules is not homogeneous in galaxies. Beside the 2.6 mm 
^^CO J = 1 ^ emission line, other more transparent molec- 
ular tracer like ^^CO J = 2 ^ 1 (Oka et al. 1998a), ^^CO and 
CS (Bally et al. 1987), as wefl as HCN (Jackson et al. 1996; 
Lee 1996), have also been used to map dense gas like in the 
Galactic centre region. 

The interpretation of features in the£ — V diagrams is sub- 
ject to many difficulties. First, since only the radial component 
of the velocity is measured, the in-plane velocity field can be re- 
covered directly only under very restrictive symmetry assump- 
tions like axisymmetry. Second, some features may not trace 
real spiral arms but be artifact of velocity crowding along the 
line of sight (Burton 1973; Mulder & Liem 1986). Third, fea- 
tures elongated in the velocity direction can be understood ei- 
ther by real structures elongated along the line of sight or by 
spatially localised emission with a very large velocity spread, 
as would result for example in a supernova explosion or in vi- 
olent shocks. Finally, some regions in the £ — V space overlay 
the emission from several distinct sources, in particular the low 
velocity emission from the Galactic centre region is partly hid- 
den by the emission from the surrounding disc spiral arms. This 
problem can however be addressed considering the latitudinal 
gas distribution or resorting to very dense gas tracers. Optically 
thick regions are also strongly affected by absorption. 



3. Numerical methods 

The initial models are evolved using a composite A^-body and 
hydro code developed by the Geneva Observatory galactic dy- 
namic group, originally provided by D. Friedli but considerably 
modified in the present application. At each time step, the grav- 
itational forces on all particles (gas included) are derived by 




Fig. 3. Example of a double grid with Nr = 10, A^^ = 12, 
Nz = 13 and = 3. The solid points and dotted lines are 
the points and meshes of grid A, and the circles and solid lines 
those of grid B. 



the particle-double mesh method detailed in Sect. 3. 1 , the pres 



sure and viscous forces on the gas particles by the Lagrangian 
smoothed particle hydrodynamics (SPH) technique described 
in Sect. and finally the phase space coordinates of each 
particle are updated by integration of the equations of motion 
according to the algorithm presented in Sect. 

3.1. Gravitation with PM'^ 

To compute the gravitational forces on the particles, we have 
resorted as in paper I to the particle-mesh (PM) technique with 
polar-cylindrical grid geometry and variable homogeneous el- 
lipsoidal kernel for the softening of the short range forces 
(Pfenniger & Friedli 1993). However, instead of using a single 
grid with constant vertical resolution, inappropriate to model 
efficiently the high density contrast between a thin disc and 
an extended dark halo, we have introduced a double embedded 
grid structure. The higher resolution grid (A) has Nr x A/'^ x 
cells, with a maximum radial extent i^max and a vertical spac- 
ing Hz, i.e. z^^^ = {Nz — l)/2 • Hz. The lower resolution 
grid (B) has half the number of cells in each dimension of the 
plane, skipping every second R- and ^-grid point of the finer 
grid, and the same number of cells vertically but with a spacing 
of Mz • Hz, i.e. z^^^ = Mz • ^max- Hence grid B has the same 
radial extent than grid A, but is larger in the vertical dimension. 



R. Fux: 3D self-consistent A^-body barred models of the Milky Way 



5 



Both grids have constant vertical spacing, thus allowing to keep 
the advantage of the Fast Fourier Transform algorithm in this 
dimension. Here Nz designates the number of active cells in z, 
corresponding to half the value required by the doubling-up 
method (Hockney & Eastwood 1981). The vertical grid param- 
eters are chosen to satisfy: 

Nz = l^2k-Mz, A: integer > 2, (1) 

ensuring that (i) there are grid points of each grid in the plane 
2: = 0, (ii) the horizontal planes z = ^z^^^ delimiting verti- 
cally grid A also contain points of grid B, and (iii) grid B has at 
least four cells within grid A in the z dimension. Figure ^ gives 
an example of double grid fulfilling all these conditions. 

In this particle-double mesh (or PM^) method, the mass as- 
signment to the grid cells is repeated twice (once for each grid), 
and the total potential involves the evaluation of three sub- 
potentials: <I>A and <I>B generated by the mass inside grids A 
and B respectively, as well as <l>c generated by the mass within 
grid B but outside grid A. The potential is computed on 
grid A, at high resolution, and the other potentials on grid B. 
The total potential then amounts to <I>a + in the region of 
grid A, which requires an interpolation of to grid A, and 
simply to <I>B outside this region. The increase of CPU time 
owing to the triple potential evaluation is amply compensated 
by the reduction by a factor ~ of the effective number of 
cells relative to a single grid with the same resolution as grid A 
and the same size as grid B. The transition planes between the 
two grids are associated with numerical discontinuities in the 
force components, although quite small in practice (see Fig. ^. 
To smooth them out, the forces in the first and last cells of 
grid B within grid A are taken as a linear combination so that 
the forces at z = :^z^^^ coincide with the forces of grid B and 
those at z = ±(^max ~ MzHz) with the forces of grid A. 

The grid parameters adopted in our simulations are given in 
Table! and the resulting gravitational resolution achieved with 
the high resolution double grid is depicted in Fig. ||. 

3.2. Gas hydrodynamics with SPH 

The pressure and viscous forces acting on the gas particles are 
derived by three-dimensional SPH. We will not repeat here the 
principle of this very popular Lagrangian method for solving 
Ruler's equation of motion (see Benz 1990; Monaghan 1992 
and Steinmetz 1996 for reviews), but just give the necessary 
specifications about this part of our code. 

The relevant form of the Euler equation writes: 

dv dv , , VP „ 

-j7 = ^^i^-^)^ = n - (2) 

dt at p 

where v{r) is the velocity field, p{r) the spatial mass density 
and P{r) the pressure of the gas. H is the viscosity term and 
— the gravitational force as obtained by the PM^ method. 

The internal properties of a classical fluid are fully de- 
scribed by three macroscopic independent variables of state. 
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Fig. 4. Typical z-behaviour of model dynamical properties cal- 
culated with the high resolution double grid, i.e. Nr = 62, 
= 64, Nz = 121, Mz = 6, i^max = 50 kpc and 
Hz = 60 pc (the mass distribution and the grid refer to model 
11012000). From top to bottom and in linear scales: mass den- 
sity (p), potential ($) and the three force components (Fr, F^f) 
and Fz). The full lines are the results from grid A (based on 
$A + ^c). and the dashed lines those from grid B (based on 
$b)- The vertical dotted lines indicate the limits z = ±^max 
of grid A. The forces considerably differ from one grid to the 
other near the plane z = 0, but almost confound in the transi- 
tion zone. 

like for example the density and the pressure appearing explic- 
itly in Eq. (||), and the internal energy per unit mass u{r). The 
density is evaluated directly at each particle position Vi via: 

Pi p{ri) ^mjW{ri - rj,h), (3) 

where A^g is the number of gas particles, rrij the individual 
mass of the particles, W{r^h) the kernel function and h the 
smoothing length. In SPH, this equation is equivalent to the 
continuity equation expressing the conservation of total mass. 

A first relation between the variables of state is given by the 
equation of state assuming a perfect gas, i.e. P = (7 — l)p • i^. 
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where 7 is the adiabatic index of the gas, which from statistical 
mechanics is given by 7 = 1 + 2//, with / being the number 
of freedom degrees of the gas constituents (/ = 3 for atomic 
hydrogen and / = 5 for molecular hydrogen). A second re- 
lation comes from the energy conservation equation, which in 
its general form involves cooling and heating functions. Since 
these functions are poorly constrained by theory, we have sim- 
ply resorted to the isothermal assumption, i.e. the internal en- 
ergy of the gas is at any time and everywhere the same. This 
assumption, partly justified by the homogeneous velocity dis- 
persion of the warm gas component in external galaxies (e.g. 
van der Kruit & Shostak 1984), means that cooling and heat- 
ing exactly cancel each other and that the energy released by 
the shocks is instantaneously radiated away. Hence, introduc- 
ing the (constant) sound speed Cg = [7(7 — 1) • ix]^/^ to sub- 
stitute the internal energy, the equation of state for particle i 
reduces to: 



P. 



7 



• Pi, 



(4) 



which is degenerated in sound speed and adiabatic index, 
depending only on the ratio Cg/7. With this last equation, 
Benz's (1990) SPH version of the pressure term in Euler's 
equation simply becomes: 



—) ^--■y,rnj(-^-)v,W{r,-rj,h). (5) 



The empirical artificial viscosity is calculated exactly as in 
Benz (1990). It consists of a superposition of the "bulk" and 
the von Neumann-Richtmyer viscosities, with the standard pa- 
rameters set to a = 1.0 and (3 = 2.5, and takes into account 
Balsara's (1995) correction to avoid energy dissipation in pure 
shearing flows. 

We have used the spherical spline kernel because of its fi- 
nite spatial extension: 




< 1^ < 1, 
l<w<2, 
w>2, 



(6) 




Fig. 5. Various spatial resolutions of the high resolution simu- 
lations presented in this paper, as a function of galactocentric 
radius. The solid, short dashed and long dashed lines respec- 
tively represent the vertical, azimuthal and radial semi-axes 
of the gravitational homogeneous ellipsoidal kernel associated 
with grid A, and the points are the smoothing lengths of a ran- 
dom selection of 10% SPH particles in model 11012000. The 
gas smoothing length is close to the vertical gravitational reso- 
lution. 



35 

§ neighbours 



Fig. 6. Typical distribution of the number of neighbouring SPH 
particles, obtained by adjusting at each time step the smoothing 
lengths according to Eqs. (|^) and (^. The data refer to model 
11012000. 



where w = \r\/h. To increase the hydro resolution in high den- 
sity regions like shocks, the particles are assigned individual 
smoothing lengths hi such that their number of neighbouring 
particles Ni always remains close to a fixed number No . Two 
particles i and j are defined as mutual neighbours Hi ^ 2 ^i^d 
|ri — Tj I < 2hij, where hi^ = {hi -\- hj)/ 2 is the symmetrised 
smoothing length which should enter the kernel evaluations in 
Eqs. (||) and (||) to ensure momentum conservation. 

At each time step, the smoothing lengths are updated con- 
sidering the general three-dimensional scaling law: 



hi_ 

ho 



^No^l Pi 



1/3 



(7) 



where ho and po are constants, and +1 is added to account for 
particle i. Following Benz (1990), differentiating this formula 
and substituting the continuity equation yields: 



hi = 



— - -h 

dt ~ 3 ' 



1 dNi 
Ni + 1~^ 



(8) 



Benz did only consider the variation of hi ensuring the con- 
stancy of the neighbour numbers, and hence did not include 
the term in Ni. However, in practice, numerical fluctuations 
of the Ni's can carry these number outside a reasonable range 
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around A^o- To overcome this problem, we simply damp the 
systematic departure of the A^^'s from No setting: 



dNi 



r]At 



(9) 



in Eq. and integrate the resulting equation along with the 
equations of motion (see Sect. |33| ). The parameter rj con- 
trols the damping rate per time step (At) and should be sig- 
nificantly greater than 1 to avoid abrupt discontinuities in the 
non-gravitational forces. In all simulations we have adopted 
A^o = 35 and r] = 5. Figure || shows the typical SPH reso- 
lution achieved in evolved models and Fig. ^ an example of 
the neighbour number distribution. The standard deviation of 
this distribution usually amounts to only about 1.5 neighbours. 
The velocity divergence in Eq. (|8|), which is also needed for the 
viscous forces, is estimated as in Benz (1990). 

The initial smoothing lengths of the gas particles, confined 
near the plane z = 0, are derived from the two-dimensional 
scaling relation: 



hi{t to) 



I NoTUi 



(10) 



where ^g{R) is the surface mass density profile of the gas com- 
ponent and to its switch-on time (see Sect. ^. This method un- 
fortunately leads to a large spread in the initial neighbour num- 
bers. To adjust the hi's without evolving too much the system, 
the simulations with live gas are therefore started with a time 
step much smaller than in the subsequent stabilised regime, 
i.e. At = 0.001 Myr for the large simulations (1-series) and 
At = 0.01 Myr for the smaller ones (s-series). 

3.3. Integrator 

To integrate the equations of motion, we have applied a semi- 
adaptative second order accurate algorithm where the (Carte- 
sian) phase space coordinates, contrary to the standard leap- 
frog scheme, are evaluated at synchronous times. At each time 
step the positions Vi and velocities Vi of the particles are mod- 
ified according to (e.g. Hut et al. 1995): 
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(11) 

(12) 
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where the 's are the accelerations and At is the time step, and 
where the last equation is used to update the smoothing lengths 
of the SPH particles. The indices n and n + 1 refer to the values 
at time T and = T + Ar respectively. 

The first equation handling the positions is not exactly re- 
versible in time. However, we have repeated simulation mOO of 
paper I replacing successively the leap-frog integrator by this 
algorithm and a second order Runge-Kutta-Fehlberg algorithm 
(RKF; Fehlberg 1968) as in Friedli's original code, and found 



Leap — frog 
This paper 
RKF 
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t [Gyr] 

Fig. 7. Variation of the total energy (top) and of the total angu- 
lar momentum about the rotation axis (bottom) for simulations 
started from the same initial conditions but evolved using three 
different second order integrators with the same constant time 
step; AE{t) = E{t) - E{0) and AL^{t) = L^{t) - L^(0). 
The steeper part of the curves are associated with the epoch of 
bar formation. 



that, for the same constant time step At = 0.1 Myr, the conser- 
vation of the total energy (and of the total angular momentum) 
is similar for the present algorithm and for the leap-frog, but 
much worse for the RKF (AE/\E\ larger by a factor ~ 20, 
see Fig. Also, the bar begins to form much sooner with the 
RKF integrator, Sitt ^ 800 Myr instead of ^ 1600 Myr with 
the other two integrators, which is a sign of higher numerical 
noise amplification. 

For the gaseous particles, the pressure and viscous forces 
contributing to a/^^^ both involve velocities which are not 
known a priori, and similarly h^~^^ depends on h^~^^ (through 
Eq. ^). Therefore the non-gravitational part of a-"^^^ and 
/i^^^ are calculated with the first order predictors: 



-a-Ar, 

/i?Ar, 



(14) 
(15) 



which fortunately preserve the second order accuracy of r'/^^^, 
and h'l^K 

To temporally resolve the high density shocks in the gas 
component, we use an adaptative time step procedure inspired 
from Friedli (1992). After step n, the next time step At"^+^ 
is estimated from the maximum relative contribution per unit 
time of the second order terms to the integrated quantities: 



(16) 
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where 



.n+l\2 



, n+1 



' lim 



1 K 

2 



h?\ 



(17) 
(18) 
(19) 



and where rum = hum and vnm = ^lim/At"^ are destinated 
to prevent divergences and a drastic drop of the time step. The 
value of h\im is fixed to 0.01 pc in the large simulations and to 
0.02 pc in the smaller ones, corresponding approximatively to 
the smallest dimension of the PM^ cells. Ideally, to guarantee 
a relative uncertainty of the same order or smaller than a given 
tolerance Etoi, it suffice to take At"^+^ = £^toi/emtx- However, 
a more advisable prescription is: 



Ar+^ = MaxjMin • A^, Ato] , At^in} , 



with 
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(20) 

(21) 
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and Qfmax = 2. The square root over (3^^^ softens the vari- 
ations of the time step and amax puts an upper limit on its 
growth. Moreover the time step is maintained in the range 
Atmin < At < Ato, with Atmin = 10"^ Myr and Ato = 
0.1 Myr. The upper boundary is a fraction of the time needed 
by the stars to cross the central high resolution PM^ cell in the 
steep nuclear potential. This time- stepping criterion is applied 
only to the gas, i.e. the loop on z in Eq. ( [l^ ) is restricted to the 
gas particles, and ej^+^ general controlled by the smooth- 
ing lengths. 

\f [3 < 2/3 and At^ > Atmin, the tolerance is considered 
as not respected in the last integration step. In this case, the old 
values of the integrated quantities are restored and integrated 
again with a time step Af = Max{ V^I^^ • Af, Atmin}. 
This rejection procedure unfortunately requires to store the old 
positions of all particles. 

In the simulations with live gas, the tolerance has been set 
to Etoi = 10~^, and in the simulations with fixed gas compo- 
nent, a constant time step At = Ato has been applied. 

4. Initial conditions 

The radial distribution of molecular gas in the Milky Way, as 
traced by CO, is essentially confined inside the solar circle, 
with a marked hole between 1 and 3 kpc and a strong central 
concentration associated to the nuclear ring/disc, whereas the 
HI presents a rather constant surface density extending far be- 
yond Ro and an abrupt decline towards the centre below 3 kpc. 
The total gas mass within Ro is estimated to 4.2 x 10^ M©, 



comprising 2.3 x 10^ Mq of H2 (with a large uncertainty 
owing to the poorly known CO to H2 conversion factor) and 
0.7 X 10^ Mq of hi (Scoville 1992), plus 28% of helium and 
metals. For the gaseous disc, we have used the following initial 
mass distribution: 



Pg{R,z) 



{27lf/^(7^pR 



exp 



R^ 



(23) 



with cr = 5 kpc, p = 0.017 and a total mass Mg = 5 x 10^ M©, 
resulting in 3.6 x 10^ Mq within the solar circle. Both radial 
and vertical profiles are Gaussian, and the disc is linearly flar- 
ing with radius, as observed in HI from a few kpc of the centre 
out to at least 2Ro (Merrifield 1992), achieving a thickness of 
136 pc at = Ro. The fast radial decline is aimed to spare 
SPH particles in the outer regions and hence increase the spa- 
tial resolution near the centre at fixed number of particles. The 
observed gas deficit between 1 and 3 kpc needs not to be repro- 
duced since shocks will naturally deplete this region. The gas 
particles initially have pure circular velocities derived from the 
axisymmetric part of the total potential at z = 0. 

The stellar and dark mass is divided into the same three 
components as in paper I: a stellar nucleus- spheroid (NS), a 
stellar disc and a dark halo (DH), with their initial axisymmet- 
ric mass distribution described by the same analytical formulae, 
except for the disc. In paper I we indeed noted that the adopted 
double exponential discs evolve into discs with too less mass 
outside the bar region according to the COBE/DIRBE near-IR 
and bulge microlensing data, and possibly an excess of mass 
in the central region when compared to the HI terminal veloc- 
ity constraints. Instead of increasing the disc scale length, we 
choose here to soften the initial central mass density taking: 



pd{R,z) = g{z/hz) 



27rhj^hz [exp (-1) + exp (-1/2)] 
exp[-^(l + ^)] if R<hR, 
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li R>h 



(24) 



where Md is the total disc mass, Hr the scale length, hz the 
scale height and g{Q the normalised vertical profile. The sur- 
face density remains exponential in the external disc, but con- 
tinuously and differentiably joins an inner Gaussian distribu- 
tion at R = Hr, with a central value reduced by 40% relative 
to the purely exponential case. To compensate for the enhanced 
spatial density near the plane z = caused by the additional 
gas component, we also replace the exponential vertical profile 
of paper I by van der Kruit's (1988) profile: 



g{C) (X sech^/" 



-c 

2^J 



(25) 



with n = 7, between the exponential (n = 00) and isothermal 
(n = 1) cases. 

The choice of the parameters are based on the simulations 
performed in paper I, giving a strong weight to the xi orbits 
versus HI terminal velocities test. The best models regarding 
this test are m06t4600 and m04t3000, whose initial axisym- 
metric models share the following interdependent properties: 



R. Fux: 3D self-consistent A^-body barred models of the Milky Way 



9 



Table 1. List and characteristics of the simulations presented in this paper. Nr, N^f^, Nz and Hz are the grid parameters as defined 
in Sect. 3.1. Both low and high resolution double grids have Mz = 6, i^max = 50 kpc, z^^^^ = 3.6 kpc and z^^^^ = 21.6 kpc. # 



is the number of particles, Cg the sound speed of the gas, to the gas switch-on time and tend the time to which the simulation is 
performed. The points refer to the former lines. 



Simulation 


Nr 






[pc] 


#NS 


# disc 


#DH 


# gas 


to [Myr] 


tend [Myr] 


Cs [kms ^] 


sxx 


32 


32 


73 


100 


83 444 


178 022 


225 284 







5000 




slO 
















20000 


3200 




10 


Ixx 


62 


64 


121 


60 


625 828 


1 335 167 


1 689 629 











110 
















150000 


1800 


2275 


10 


120 






















20 


110' 


















2400 


2775 
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(i) a ratio of disc over NS mass of 0.8 within the spheroidal 
volume s <3 kpc (where 5^ = i?^ + z^/e^ and e = 0.5), (ii) a 
total NS+disc mass of 2.4 x 10^^ Mq in the same volume and 
(iii) a circular velocity of 190 kms~^ just after the very steep 
central rise of the rotation curve, taking into account the veloc- 
ity scales adjusted to the observed stellar velocity dispersion 
in Baade's window (see Table 3 of paper I). The DH param- 
eters, a and hz are fixed as in the reference simulation mOO, 
and the other parameters are adjusted to the former constraints 
and to a disc surface density on the solar circle of 60 Mq pc~^, 
yielding Mns = 2.48 x 10^° M©, Md = 4.6 x 10^° M© and 
hn = 3.2 kpc. The mass density is softly truncated as in pa- 
per I, at radius Rc = 38 kpc and over a width S = b kpc, and 
the initial kinematics rests on the same relations between the 
velocity moments as in simulations mOO-mlO, except that the 
asymptotic velocity anisotropy Poo =0.68 and the transition 
radius To = 5.3 kpc for the NS component. 

5. Time evolution 

Table |l| gives an overview of the simulations described in this 
paper. Two series of simulations have been realised: the first 
one with moderate number of particles and spatial resolution 
(series "s", small simulations), and the other one pushing these 
quantities such as to exploit about half of the memory resources 
of the most powerful computers available at Geneva University 
(series "1", large simulations). In each series, we first integrated 
the initial axisymmetric model for 5 Gyr keeping the gaseous 
component fixed (simulations "xx"), and then we relaxed the 
gas particles at different intermediate times to . The gas is not 
evolved from the beginning because the non-inclusion of star 
formation produces an excessive accumulation of gas in the 
central region due to the torques of the non-axisymmetric grav- 
itational potential, raising the rotation curve and leading to 
a premature destruction of the bar (Friedli & Benz 1993) or 
even preventing its formation. Instead of introducing an ar- 
tificial gas recycling procedure, the time consuming gas-live 
simulations of the 1- series were integrated only over a few 
100 Myr. To damp the initial disequilibrium of the gas owing 
to its circular kinematics in the already barred potential, the 
non-axisymmetric part of the potential is progressively and lin- 



early brought to its nominal value in 75 Myr, i.e. roughly half 
a rotation period of the bar. 

In the 1- series, the gas has been released at two different 
times after the formation of the bar, at to = 1.8 and 2.4 Gyr. 
In the first case, two values of the sound speed have been ex- 
plored, Cs = 10 and 20 kms~^ (simulations 110 and 120 re- 
spectively), and in the second case only Cg = 10 kms~^ has 
been retained (simulation 110'). In the s-series, many runs with 
live gas have been performed, releasing the gas every 400 Myr 
and each time with two different sound speeds, but only the 
one mentioned in Table |l]will be discussed here. The adiabatic 
index of the gas is set to that of neutral hydrogen, i.e. 7 = 5/3. 
Taking the vertical gravitational resolution Hz as a lower phys- 
ical limit for the SPH smoothing length, the maximum gaseous 
density that can be reasonably modelised in the 1- simulations is 
- W{0,Hz) • No/Ng - Mg ^ 2 MqPc-^ (see the previous 
sections for the meaning of the symbols), which is at least one 
order of magnitude below the density of the Galactic molecular 
gas, but sufficient to describe the warmer neutral phase with a 
typical sound speed of 10 kms~^. However, Cowie (1980) has 
argued that a system of molecular clouds may be treated as a 
classical fluid with a sound speed equal to the mass averaged 
cloud velocity dispersion, which also amounts to 10-20 km s~^ 
in the Milky Way. Hence our simulations are expected to dis- 
play properties of both the HI and H2 medium. 

The SPH particles all have the same time-independent 
mass. The number of particles in each luminous component 
is always such that the mass per particle is the same as for the 
gas component to minimise relaxation effects, and the num- 
ber of DH particles corresponds to a mass per particle three 
times larger than for the other components. The typical num- 
ber of SPH particles within R = 8 kpc (in initial units) and the 
corotation circle in the high resolution simulations is 10^ and 
5.5 X 10^ respectively. 

The simulations are run in a completely self-consistent 
way, without imposing any symmetry and taking into account 
all gravitational interactions between the mass components. In 
particular, the gas feels his own gravity and interacts with the 
stellar spiral arms. Moreover the bar parameters are not arbi- 
trarily chosen but automatically and naturally adjust according 
to realistic dynamical constraints. The calculations have been 
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Fig. 8. Time evolution of the face-on disc+NS surface density in simulations Ixx and sxx with fixed gas component. The distances 
are in kpc (initial units) and the density contours are spaced by a constant interval of 0.5 magnitude. The cross in the Ixx frames 
indicates the position of the centre of mass, which coincides with the origin of the coordinates system. 
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Fig. 9. Radial, x- and ^-displacements of the stellar density 
centre (dc) with respect to the centre of mass in simulation Ixx 
(thin line) and in the live gas simulations 110 and 110' (thick 
lines). 

done on Sparc Ultra and Silicon Graphics computers, each with 
over 1 Gbyte central memory. 

5.7. Stars 

Figure || shows the whole face-on evolution of the simula- 
tions Ixx and sxx with rigid gas component. Although the ini- 
tial conditions of these simulations are drawn from exactly the 
same phase- space density function and the average number of 
particles per grid cell is about the same, the evolution clearly 
depends on the adopted resolution and number of particles: in 
simulation Ixx the bar forms much more rapidly than in simu- 
lation sxx, around t = 1.2 Gyr instead of 3.2 Gyr in the latter, 
and is of larger extent. Moreover, the contours of the bar be- 
come rounder close to the centre, as observed in many external 
barred galaxies (e.g. Ml 00 in the near-IR), and the surrounding 
disc tends to a more flattened radial profile. An explanation for 
the delay of the bar formation could be that higher resolution 
simulations can catch smaller density fluctuations and hence 
favour the growth rate of asymmetries. It is not clear whether 
a convergence of properties with increasing resolution has been 
achieved in simulation Ixx. 

Another relevant dynamical aspect distinguishing the large 
simulations from those of the s-series is the offcentring of the 
stellar bar (Fig. ^. At t 800 Myr, the density centre starts to 
deviate from the global centre of mass and wanders around it. 
The maximum amplitude of the displacement reaches ~ 800 pc 
at t = 2.9 Gyr, and the revolution frequency of the density cen- 
tre amounts to 20 — 30 km s~^ kpc~^ . Offcentred bars are com- 
monly observed in external galaxies (see Colin & Athanassoula 
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Fig. 10. Pattern speed ftp {t) of the bar in simulation Ixx with 
fixed gas (thin line) and in simulations 110 and 110' with live 
gas (thick lines), in initial units. The pattern speed is derived as 
in paper I by diagonalisation of the momentum of inertia ten- 
sor, except that the latter is computed relative to the offcentred 
density centre instead of the centre of mass. 

1989; Block et al. 1994; other references in Levine & Sparke 
1998) and reported in numerical simulations of galactic discs 
(e.g. Miller & Smith 1992). At least half of all spiral galaxies 
have lopsided light distribution (Schoenmakers 1999; see also 
Rudnick & Rix 1998). Figure | confirms Weinberg's (1994) 
conclusion on the persistent nature of the phenomenon. How- 
ever, it might be that the polar grids used to compute the grav- 
itational forces amplify the density centre offcentring, as they 
artificially produce an exponential instability of the position of 
the centre of mass (Pfenniger & Friedli 1993). The same sim- 
ulation repeated with a three-dimensional Cartesian grid could 
result in a lower amplitude oscillation. 

The pattern speed of the bar in simulation Ixx decreases 
roughly exponentially from 50 km s~^ kpc~^ at t = 1.2 Gyr to 
30 kms~^ kpc~^ at t = 5 Gyr (see Fig. |lo|). The face-on axis 
ratio b/a of the bar is about 0.6, i.e. close to the upper limit 
derived from the lower resolution simulations of paper I. 

The phase space coordinates of the stellar particles have 
been extracted from the simulations every 25 Myr. Adjusting 
the location of the observer (Sun) using the COBE/DIRBE 
dust subtracted K-band map as in paper I is complicated by 
the density centre offcentring: in addition to the bar inclina- 
tion angle, the relative galactocentric distance of the observer 
and the mass-to-light ratio, further parameters are needed for 
the direction of the Galactic centre in the models. Here we 
have simply applied the standard method of paper I to mod- 
els with weakly offcentred bars, assuming that the Galactic 
centre lies at the centre of mass and taking a finer Carte- 
sian grid with A£ = Ab = 0.5° to compare the data and 
model fluxes. In the notations of paper I (see also Sect. ^ and 
for A/'pix = 1200, the best fit location parameters for model 
lxxtl950 are = 9.7 and (fo = 32°, with a mean quadratic 
relative residual 7^^ = 0.539%. A comparison of this model to 
the COBE data is shown in Fig. |ll| 
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0.5 1 1.5 2 2.5 Logio(M0pc ^) 

Fig. 12. Face-on view of the gas flow evolution in simulations 110 and 120, which differ only by the value of the sound speed Cg. 
Each frame is 20 kpc on a side in initial units. The dotted lines indicate the stellar surface mass density contours spaced by 0.75 
magnitude. 
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Fig. 11. Result of an adjustment of model lxxtl950 to the dust 
subtracted COBE K-band image by the same technique as de- 
scribed in paper I. The solid lines show the observed contours 
spaced by half a magnitude and the dashed lines the corre- 
sponding model contours, assuming a constant mass-to-light 
ratio. The correction for extinction fails below the horizontal 
line. 



5.2. Gas 

The time evolution of the gas flow in simulations 110, 120 and 
110' is illustrated in Figs. |l^ and |l3|. In each simulation the 
gas flow rapidly becomes non-axisymmetric, forming transient 
spiral arms, shock fronts and a nuclear ring of X2 orbits accu- 
mulated near the inner Lindblad resonance. Two kinds of spiral 
structure can a priori be distinguished in the bar region (see for 
instance frame t = 2075 Myr of simulation 110 in Fig. [l^ or 
frame t = 2775 Myr in Fig. [l|): 

- The axis shocks, or off-axis shocks, which lead more or less 
the bar major axis and join the nuclear ring. These shocks, 
also appearing in many other hydrodynamical simulations 
(e.g. Athanassoula 1992), can be roughly understood on the 
basis of the xi closed orbit family in the rotating frame of 
the bar (Binney et al. 1991; Morris & Serabyn 1996), under 
the approximation of a rigid potential. Far from the centre, 
the gas moves along this main orbit family because the vis- 
cous forces dissipate any libration energy around periodic 
orbits. The same forces also cause the gas to switch pro- 
gressively to ever lower energy orbits and thus to approach 
the centre. Below a critical value of the Hamiltonian, the 
xi orbits develop loops at their apocentre where the gas dis- 
sipates some of its streaming energy by collision. The gas 
then leaves these periodic orbits to follow more radial non- 
periodic orbits passing round the nuclear ring and striking 
the gas falling symmetrically from the other side of the bar. 
The axis shocks result from the velocity difference between 
the two streams, when it exceeds the sound speed. Part of 
the falling gas may also collide with the central X2 orbits 
and be directly absorbed by the nuclear ring. Axis shocks 
have been detected in the velocity field of several external 




Fig. 13. Face-on view of the gas flow evolution in simulation 
110'. Each frame is 20 kpc on a side in initial units. The gray 
scale and the dotted contours are as in Fig p2l 



barred galaxies (see Sect. 6.3) and seem to be associated 
with the prominent dustlanes leading the bar in these galax- 
ies. The dust grains are strongly concentrated behind the 
shock fronts and thus produce the typical extinction signa- 
tures. Some prototypical galaxies with offset dustlanes are 
NGC 1300, NGC 1433, NGC 1512, NGC 1530, NGC 1365 
and NGC 695 1 , which are all SBb or SBbc type galaxies 
(see e.g. Sandage & Bedke 1988), as the Milky Way. The 
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self-gravity of the gas certainly plays an important role in 
holding together the shocked gas. 
- The lateral arms, which roughly link the bar ends avoiding 
the nuclear ring by a large bow, and generally correspond 
to the inner prolongation of spiral arms in the disc. The 
gas moves almost parallel to these arms and finally meets 
the axis shocks. According to Mulder & Liem (1986), the 
3-kpc arm is almost certainly of this kind. When a lateral 
arm forms, the part close to the axis shock where the gas 
runs into is located well inside the bar. Then the arm moves 
outwards along the axis shock and progressively dissolves 
as it approaches corotation (this will be illustrated more 
quantitatively in Sect. 0). Such dissolution may be partly 
linked with the decreasing gravitational resolution with ra- 
dius in the simulations. The gaseous lateral arms in our 
simulations resemble the innermost, sometimes ring-like, 
stellar arms in external barred galaxies, e.g. NGC 1433, 
NGC 4593, NGC 6951, NGC 3485, NGC 5921, NGC 7421, 
which are also SBb or SBbc type galaxies (Sandage & 
Bedke 1988). 



lateral arms pass away from it. 

Increasing the sound speed, i.e. the pressure forces relative 
to the gravitational forces, yields as expected a smoother gas 
distribution and softens the sharp corners in the spiral arms (see 
Fig. ^2|). The nuclear ring slightly shrinks and the axis shocks 
occur closer to the major axis of the bar, in agreement with the 
SPH results of Englmaier & Gerhard (1997). Note that some 
sub- structures also seem to develop within the nuclear ring for 
Cs 20kms-^ 

The gas flow never reaches a stationary state and is most 
of the time asymmetric, as is often observed in external disc 
galaxies. In particular, the lateral arms are rarely symmetric 
both in position and surface density, and the axis shocks consid- 
erably vary in shape. Sometimes, only one lateral arm or axis 
shock can be distinguished, and sometimes these structures are 
doubled, appearing twice on the same side of the bar (see e.g. 
11012175 in Fig. |l2|for double lateral arms). It may even hap- 
pen, as is the case in frame t = 1975 Myr of simulation 110 
(Fig. [l^), that the axis shocks receive a sufficient impulse from 
the lateral arms to transform themselves into such arms. The 
nuclear ring, owing to its gravitational coupling to the stellar 
components, intimately follows the bar density centre in its os- 
cillations around the centre of mass. The inclusion of gas pre- 
vents the pattern speed of the bar to slow down (Fig. |0f see 
also Friedli & Benz 1993; Berentzen et al. 1998): after a sub- 
stantial initial readjustment, ftp remains very nearly constant 
over the short duration of the live gas simulations. 

In the s- series, the moderate bar considerably weakens once 
the gas is switched-on. However, in simulation slO it survives 
for a sufficiently long time to reveal an interesting asymmetric 
process in the gas flow (Fig. |lj): the gas seems to rarefy alter- 
natively on each side of the bar intermediate axis according to 
a cycle of roughly 125 Myr. As a large amount of gas flows on 
one side, the other appears almost devoided of gas, and when 
the gas of the former side reaches the axis shock, the situation 
reverses and the previously depleted side is filled. During this 
cycle, the gas density on the bar flanks can vary from single to 
triple. 

The non-stationarity of the gas flow in the 1- simulations is 
not a simple unachieved readjustment due to non-equilibrium 
in the initial conditions and the rather short integration time, 
since the smaller simulations with live gas, integrated much 
longer in time, themselves always remain time-dependent. Out- 
side the bar region, the stellar spiral arms have a differential 
pattern speed, in general lower than the bar (see e.g. Sellwood 
& Sparke 1988), and drive the gaseous arms which repeatedly 
wound and dissolve. Many SPH simulations carried in rigid 
barred potentials seem to finally settle in a stationary spiral 
structure rotating at the same angular speed as the bar. Such 
results are an artificial consequence of the imposed fixed back- 
ground potential. 



This classification should only be considered as a first order 
guide. The main distinction between the two structures is that 
the axis shocks intersect or brush the nuclear ring, while the 



6. Interpretation of the i — V features 



Model £ — V diagrams depend on the location of the observer 
in the galactic plane, which is described by his galactocentric 
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Fig. 14. Asymmetric alternated cycle of the gas flow on the bar flanks in the small simulation slO. The full and dotted lines 
indicate the gas mass inside each of the two regions corotating with the bar, located within 3 kpc from the centre and at least 
500 pc away from the bar major axis. The cycle is especially marked in the last 700 Myr of the simulation (in which the bar never 
gets significantly off centred). 



Table 2. Reference points overlayed in the £ — V diagram 
movies to locate the main observed features indicated in Fig. |l] 
and the Carina arm. 



Reference 




V [kms-^] 


Positive bright velocity peak 


3 


270 


Negative bright velocity peak 


-2 


-220 


Intersection of the 3-kpc arm with the ^-axis 


12 





Knee of the 3-kpc arm 


-17 


-125 


Knee of the 135-km s~^ arm (Clump 1) 


-5 


100 


Maximum emission of the connecting arm 


7.5 


150 


Tangent points of the molecular ring branches 


25 


120 




30 


115 


Tangent point of the Carina arm 


-78 






Two models from these movies, 110t2066 and 110't2540, 
illustrated and confronted to the observations in Fig. [l^, have 
especially retained our attention. In both cases, the inclination 
angle of the bar is 25° and no velocity rescaling has been ap- 
plied. The circular velocity of the observer is set to the mean 
azimuthal velocity of his surrounding SPH particles, i.e. to 
203 kms-^ for 110t2066 and 197 kms"^ for 110't2540. The 
model 110't2540 is selected as one of the models which best 
reproduce the overall l — V observations within the solar circle. 
Even if the agreement is rather qualitative, this model offers a 
solid basis to interpret the data. Figure |l^ highlights the spiral 
arms which trace the ^ — F features in the model and Fig. [l^ 
shows the velocity field of the gas. Figure |l8| also provides a 
qualitative key-map to spatially locate observed i — V sources 
in the Galactic plane. 



distance Ro in initial units and the angle Lpo between the line 
joining himself to the centre and the major axis of the bar, and 
on his velocity Vo- We assume here that the observer lies in the 
plane z = and that is purely azimuthal. 

The phase space coordinates of the SPH particles have been 
stored every 2 Myr in order to realise a posteriori ^ — V dia- 
gram movies for any arbitrary values of these parameters with 
a high temporal resolution. To compare these movies with the 
observations and decide for an optimum model, some refer- 
ence points have been overlayed on the movies to localise the 
features seen in the CO and HI data (Table In each movie, 
the view point is at rest relative to the rotating frame of the 
bar, i.e. and Lpo remain constant. The contribution of each 
SPH particle to the ^ — F frames is weighted by its inverse 
squared distance relative to the observer to mimic the flux de- 
cline of point sources (all model ^ — V diagrams in this paper 
are computed this way, except the one in Fig. A direct con- 
sequence of the asymmetries developing in our simulations is 
that model l — V diagrams computed for diametrically opposite 
view points always considerably differ. 



6.1. Connecting arm 

The connecting arm is clearly identified with the axis shock in 
the near part of the bar. More precisely, this arm is build up 
by the gas clouds which have recently crossed the shock front 
at various galactocentric distances and which now collectively 
plunge towards the nuclear ring/disc with velocities roughly 
parallel to the shock front. In other words, the connecting arm 
traces the near- side branch of the Milky Way's dustlanes. Other 
models, like 110t2066 in Fig. |l5[ exhibit axis shocks with l — V 
traces resembling much more the real connecting arm than in 
model 110't2540, thus reinforcing our interpretation. 

The presence of the connecting arm feature in the observed 
^ — V diagrams can be considered as a further evidence of the 
Galactic bar. Furthermore, its rather large domain in longitude 
is relevant of an extended bar, especially if the latter is seen 
close from its major axis (see Sect. ^. The connecting arm 
feature is a real concentration of gas in space and not an artifact 
due to velocity crowding along the line of sight, as suggested 
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by Mulder & Liem (1986)[} In our simulations, the gas does 
not always trace the full length of the axis shocks and hence 
it is not surprising that the emission from the connecting arm 
appears truncated at ^ ^ 10°. The £ — V movies sometimes 
show gas lumps deposited by lateral arms in the near-side axis 
shock and moving precisely along the connecting arm feature. 
The time for the gas to travel from the endpoint of a lateral arm 
to the nuclear ring along the shock is of order 15 — 20 Myr and 
it takes about another 10 Myr for the gas not absorbed by the 
ring to encounter the opposite axis shock. The fraction of gas 
deposited in the nuclear ring is time-dependent, but has been 
estimated to 20% in a steady gas flow model of NGC 1530 
(Regan et al. 1997). The gas mass in each branch of the axis 
shocks is ~ 4 X 10 Mq or less. 

According to Fig. |l^, the axis shock in the far- side of the 
bar is predicted as an almost vertical feature in the £ — V di- 
agram, i.e. with £ ^ constant. Such a feature is indeed visible 
in the ^^CO observations (Fig. |l]) near £ = —4°, with only a 
marginal decrease of absolute longitude towards negative ve- 
locities. The longitude confinement comes from the fact that 
the shock line is nearly parallel to the line of sight and the 
velocity extension from the fact that the velocity of the gas 
along the shock rapidly increases when approaching the nu- 
clear ring/disc. 

6.2. 3-kpc and 135-kms~^ arms 

The 3-kpc and the 135-kms~^ arms are lateral arms. The ob- 
served velocity asymmetry between these two arms happens 
because the latter passes closer to the centre than the former. 
The gas associated with the 135-kms~^ arm, moving almost 
parallel to the arm (Fig. |l7|), indeed falls deeper in the poten- 
tial well of the nucleus and therefore reaches higher forbidden 
velocities before striking the dustlane shock. 

The observed velocity asymmetry of bright emission near 
the positive and negative terminal velocity peaks could have a 
related origin: the gas from the 3-kpc arm, after crossing the 
dustlane shock further out than its counterarm, starts to fall to- 
wards the nuclear ring/disc with a higher potential energy and 
therefore will approach the latter with higher velocities, pro- 
ducing an enhanced velocity peak. Model 110' 12540 however 
does not exhibit such a velocity asymmetry. This asymmetry 
could also arise from a relative radial motion between the nu- 
clear ring/disc and the LSR, produced either by the oscillations 
of the density centre (see Sect. 5J_) or by a radial velocity com- 
ponent of the LSR with respect to the Galactic centre, possibly 
induced by the bar itself (Raboud et al. 1998), or both. 

The inner branch of the molecular ring, with tangent point 
at ^ ~ 25°, is the outer extension of the 135-kms~^ arm. 



6.3. Bania's clumps 

Observations of the gas velocity field in external barred galax- 
ies have revealed velocity changes up to 200 km s~^ across the 
bar leading dustlanes, demonstrating that these dustlanes are 
associated with very strong shocks. Such velocity fields have 
been measured for example in NGC 6221 (Pence & Blackman 
1984), NGC 1365 (Joersaeter & van Moorsel 1995; Lindblad et 
al. 1996), NGC 3095 (Weiner et al. 1993), NGC 1530 (Regan 
et al. 1997; Reynaud & Downes 1998) and NGC 7479 (Laine et 
al. 1999). Figure |l9| shows how the gas velocity field is affected 
by the near- side branch of the axis shocks in model 11012066. 




^ These authors pretend that the connecting arm intervene between 
longitudes 10° and 20°, which obviously conflicts with the sources 
(van der Kruit 1970 and Cohen 1975) they mention. 



Fig. 16. Link between the spiral arms in thcx — y plane and the 
£ — V features in model 110't2540. The spiral arms and their 
£ — V traces are depicted by the true phase space coordinates 
of the gas particles, using all particles within a narrow band of 
300 pc centred on the maximum surface density curve of each 
spiral arm. In the upper plot, closer structures relative to the 
observer have been overlay ed on the more distant ones. The 
nuclear ring is not represented. 
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X [kpc] X [kpc] 

Fig. 15. Confrontation of a selection of two models with observed gas kinematics. Top: ^^CO and Hli — V diagrams integrated 
over \b\ < 2° and \b\ < 1.25° respectively; the data are from Dame et al. (1999) for the CO, and Hartmann & Burton (1997), 
Burton & Liszt (1978) and Kerr et al. (1986) for the HI. Middle: synthetic i-V diagrams of models 11012066 and 110't2540 
for a bar inclination angle (fo = 25°, including all particles within \b\ < 2°. Bottom: face-on projections of the gas spatial 
distribution in these models, rescaled such as to put the observer at (x, y) = (0, —8) kpc (0 symbol). In these units, corotation 
lies at i?L = 4.5 kpc (11012066) and 4.4 kpc (110't2540). The model on the left reproduces almost perfectly the connecting arm, 
while the model on the right provides a fair global qualitative agreement to the data. 
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The shock occurs near x = —0.55 kpc and is followed down- 
stream by a huge gas concentration. Both velocity components, 
parallel and perpendicular to the shock front, undergo an abrupt 
velocity gradient across the shock layer. The velocity change is 
however larger for the parallel component, reflecting the shear- 
ing nature of the axis shocks, and is comparable to the val- 
ues observed in external galaxies. In NGC 1530, Reynaud & 
Downes (1998) found that the velocity change increases to- 
wards the nuclear ring, but this property is hard to infer from 
our models because the shock fronts are not well resolved in 
the low density part. 

Bania's (1977) clump 2 and an other vertical feature near 
£ ^ 5.5° (see Fig. [l|) could represent gas lumps which are just 
about to cross the near-side dustlane shock somewhere between 
the 3 -kpc arm and the nuclear ring/disc: the upstream part still 
moves with the small quasi- apocentric velocities of the pre- 
shock orbits, while the downstream part has been accelerated 
to the high inwards post-shock velocities, giving rise to a steep 
radial velocity gradient (Fig. [T^) and a velocity stretch of over 
100 km s~^ in the observations. Contrary to the £ — V trace of 
the far- side dustlane shock, these features might be really con- 
centrated in space and not result from an accumulation of gas 
along the line of sight. For the £ = 5.5° lump, this interpre- 
tation is supported by the fact that all emission from the lump 
originates at nearly constant latitude, as expected from a spa- 
tially confined source, and that the part of the connecting arm 
at the same longitude as the lump appears at almost the same 
latitude as the lump itself (see Fig. g). Furthermore, this lump 
also has a small mass relative to the connecting arm and will 




X [kpc] 

Fig. 17. Velocity field of model 110't2540 in the rotating frame 
of the bar. The subtracted solid rotation does not take into ac- 
count the bar off centring. The grey curves indicate the location 
of the spiral arms. 
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Fig. 18. Radial velocity contours in model 110' 12540 as a func- 
tion of longitude and distance relative to the observer, whose 
location and circular velocity are as in Fig. [l^. The spacing 
of the contours is 20 km s~^ and the thick contour is the zero 
velocity curve (labels are given in kms~^). Only the particles 
within \b\ < 2° have been considered. Note the very sharp ve- 
locity gradient when crossing the axis shocks. 



therefore essentially adapt its momentum to that of this arm. 
For clump 2, with a mean latitude differing by more than 0.5° 
from that of the connecting arm (at similar £) and a total mass 
of nearly 10^ Mq (Stark & Bania 1986), our interpretation is 
more speculative. However, if this clump is indeed close to the 
apocentre of its orbit, it will enter very slowly the shock line, 
where gas is moving at very high speed (over 200 km s~^), and 
thus receive a significant impulse when integrated over time. 
Moreover, the clump complex may move on a kind of looped 
quasi-xi orbit and therefore will self-dissipate its energy if its 
size is comparable to that of the loops, whatever its mass. If the 
connecting arm indeed traces the near-side dustlane, the identi- 
fication of clump 2 with such a dustlane, as proposed by Stark 
& Bania (1986), is ruled out (unless the Milky Way has a dou- 
ble bar). But it should be noted that the axis shock assigned to 
the connecting arm in model 110' 12540 can produce 3.n £ — V 
trace resembling much more that of clump 2 if the bar inclina- 
tion angle is reduced to 15° (see Fig. ^2| ). 

Clump 1 , composed of several clouds which are not bound 
to each other (Bania et al. 1986), is the southern terminus part 
of the 135-kms~^ arm which penetrates the far-side dustlane 
shock. According to model 110' 12540, its dynamics should be 
rather subtle: a part of its gas is absorbed by the dustlane, result- 
ing in a huge velocity change like those described above, from 
V ^ 100 kms~^ to F < 0, and the other part, corresponding 
to the portion of clump 1 at £ ^ —5°, is gliding outwards along 



R. Fux: 3D self-consistent A^-body barred models of the Milky Way 



19 



but without crossing the shock front until apocentre is reached 
(see the magenta segment in Fig. |l6|). The momentum injected 
by this cloud into the axis shock gas bends the outer segment of 
the shock. Since the far- side dustlane is nearly aligned with the 
line of sight, the vertical feature near £ = —4° is in fact a super- 
position of dustlane gas moving along this line and of clouds 
with shock induced velocity gradients (see Fig. Clump 1 
must be a very perturbed and compressed region, and hence a 
potential site of star formation a priori. An indicator of "readi- 
ness" for star formation is given by the ^^CO J = 2 ^ 1 to 
J = 1 ^ ratio, tracing dense molecular clouds, which is in- 
deed particularly high for this clump (Hasegawa 1997, private 
communication). However, the strong shearing in the shock 
may prevent any star formation to proceed (e.g. Reynaud & 
Downes 1998). The bulk of this clump, owing to its large mass 
and impact velocity, may also cross the shock front without be- 
ing too much affected. 

6.4. Tilt of the gaseous disc 

A further observational argument supporting our interpretation 
of the dominant i — V features is the fact that the connect- 
ing arm and the portion of the 3-kpc arm at positive longitude, 
located in the near- side of the bar according to the models, 
have their maximum emission below the Galactic plane (i.e. 
at 6 < 0), and that the dustlane near £ = — 4° and the negative 
longitude part of the 135-kms~^ arm, located in the far- side of 
the bar, above this plane (i.e. h > 0; see Fig. ^. Hence struc- 
tures predicted to be spatially close to each other by the models 
are indeed found at similar latitude in the observations. 

Consequently, the gaseous disc within the central 1 — 2 kpc 
is tilted relative to the plane b = 0. Referring to Fig. ^ emission 
of the connecting arm at ^ ^ 5°, i.e. where b stops to decrease 
as £ increases, occurs near b = —1°. Approximating the dust- 
lanes by a straight line across the Galactic centre with an in- 
plane inclination angle of (po = 25° (see Sect. 0), the source of 
this emission is tilted by 6> = tan~^| simpo tan 6/ sin^| = 4.8° 
out of the Galactic plane, corresponding to a distance of about 
120 pc below this plane. Similarly, the far- side dustlane ends 
Siti ^ -4.5° with b ^ 0.8°, implying a tilt = 4.3° and a 
height of ~ +135 pc. A similar tilt is also apparent on the ^^CO 
longitude-latitude map of Dame et al. (1987; 1999), where the 
dustlanes contribute only very little to the total emission. Blitz 
& Spergel (1991) have inferred an apparent 5.5° central tilt of 
the bar major axis from balloon 2.4 fim observations of the 
Galactic bulge within \£\ < 12° and \b\ < 10°. However, 
such a tilt has not been confirmed by the more recent near-IR 
COBE/DIRBE maps (Weiland et al. 1994). In the dust sub- 
tracted K-band map (paper I), the latitude centroid as a func- 
tion of longitude, excluding the high extinction zone \b\ < 3°, 
only shows a significant tilt when regions beyond 10° from the 
Galactic plane are included, but then the tilt is likely an artifact 
due to the growing contribution of zodiacal light. The small tilt 
angle derived here is not generated by a position of the Sun 
above the Galactic plane, which is only of order 10 — 40 pc 
(Humphreys & Larsen 1995 and references therein), and con- 
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Fig. 19. Spatial density and velocity profiles along four slits 
perpendicular to the "connecting arm" in model 11012066. The 
upper and lower velocity curves respectively stand for the ve- 
locity components parallel (Vy) and perpendicular (Vx) to the 
shock front. The profiles are derived by standard SPH summa- 
tion, e.g. from Eq. (jsj) for p, and for z = 0. The vertical dis- 
tribution of the dense post-shock gas everywhere peaks within 
10 pc from that plane. 



cerns gas on a larger scale than the 180-pc molecular ring, for 
which an apparent tilt angle of ~ 6° has been estimated, also 
with negative latitude in the first Galactic quadrant (Uchida et 
al. 1994a; Morris & Serabyn 1996). 

Larger tilts of the inner ~ 2 kpc gaseous disc, of 20° — 30° 
and based on expanding or elliptical motion models, have been 
reported in the past (e.g. Cohen 1975; Burton & Liszt 1978; 
Liszt & Burton 1980). More recently. Burton & Liszt (1992) 
have updated and improved their tilted disc model into a flaring 
warp with a central disc coplanar to the Galactic plane. The 
warp is rectilinear in each radius and the tilt of its midplane 
is given by ^max • sin (^ — ^o), where (j) is the galactocentric 
azimuth relative to the Sun defined positive in the direction of 
Galactic rotation, (bo = 45° the azimuth of the line of nodes 
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and ^max = 13° the maximum tilt angle. This model predicts a 
tilt angle of about 4°at(/) = (/Po=25°,in agreement with the 
tilt derived for the Galactic dustlanes. 

The departure of the dustlanes from the plane b = does 
not appear to increase linearly with galactocentric distance, but 
seems to stabilise or gently decline at larger distance, rendering 
its description by a constant tilt somewhat oversimplistic. 

6.5. Central molecular zone 

This part of the ^ — F observations is probably the most com- 
plex and the most difficult to understand. Binney et al. (1991) 
have given a detailed description of the gas flow in a rotating 
barred potential in terms of the closed xi orbits, associating 
the terminal velocity curves in the observed I — V diagrams 
with the envelope of non self-intersecting such orbits and the 
parallelogram structure of the 1 80-pc molecular ring with the 
innermost orbit of this kind, called the "cusped xi" orbit, where 
shocks would transform most of the atomic gas into molecules 
and force the gas to plunge onto the more viable orbits of the 
X2 family. 

This model is confronted to several problems. The most 
frequently reported one (e.g. Kuijken 1996; Morris & Serabyn 
1996) is the asymmetry of the observed 1 80-pc parallelogram, 
and in particular the substantial ~ 140 kms~^ forbidden ve- 
locity near I = —0.8°, which cannot be fully accounted for by 
projection effects of the cusped xi orbit. Other problems are: 
(i) when viewed in the full £ — b — V data cube, the 1 80-pc 
parallelogram appears more as an assemblage of larger scale 
features rather than forming a distinct unity (see e.g. Fig. 4 in 
Morris & Serabyn 1996); for instance, its upper edge is ob- 
served to extend continuously well beyond the longitude range 
of the postulated parallelogram and joins the connecting arm 
through the positive terminal velocity peak, (ii) the terminal 
velocity peaks and the parallelogram cannot both be generated 
by the same cusped xi orbit because the formers occur well 
outside the longitude boundaries of the latter, (iii) the longitu- 
dinal edges of the parallelogram near ^ = — 1° and 1.5°, as- 
similated to the bar leading shocks in the model, should define 
rather clear £ — V features as in our models, yet these edges are 
more disordered than the other edges of the parallelogram, and 
(iv) the cusped xi orbit of the stellar dynamical models in pa- 
per I generally does not match the HI terminal velocity peaks; 
the models providing the best agreement always arise in young 
and not completely stabilised bars, whereas in older bars the 
cusped xi orbit presents velocity peaks at fairly larger absolute 
longitudes than observed. 

Figure^ shows that the dustlane shocks are indeed respon- 
sible for the positive and negative terminal velocity peaks, but 
that they do not coincide with the leading edges of the cusped 
xi orbit, which produces velocity peaks at higher absolute lon- 
gitude and with lower velocity amplitude. Gas on the shocks 
move along non-periodic orbits with much smaller pericentre 
than the cusped xi orbit. Thus the Milky Way's cusped xi or- 
bit is probably much larger than in the Binney et al. model, ex- 
plaining why these authors find a very small corotation radius 
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Fig. 20. Typical path in the x — y and i — V planes of an SPH 
particle in simulation 110. The points show the positions of the 
particle from t = 1994 to 2200 Myr with a constant time inter- 
val of 2 Myr. The x — y plot is in the bar rotating frame with 
the origin at the density centre and with the bar inclined by 25° 
relative to x = (like in Fig. 15 ), and the ^ — F plot is viewed 
from (x, y) = (0, —8) kpc. The particle is chosen to pass 
through the intersection of the near- side lateral arm and axis 
shock ait = 2066 Myr, i.e. through (x, y) = (-1.4, -2.4) kpc 
in the lower left frame of Fig. [l^. Note the larger velocity gap 
between the points when the particle crosses the far-side axis 
shock. 



of 2.4 kpc. The asymmetries in our models make the cusped xi 
orbit slightly uncertain. In particular, there is a small range of 
Hamiltonians where the xi orbits only loop at one extremity. 
However, studies of Galactic xi orbits generally represent the 
true Milky Way's potential by bisymmetric models and there- 
fore are also concerned with similar uncertainties. 

The gas on the dustlanes falling onto the nuclear ring from 
large distances does not necessarily merge with the ring, but 
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Fig. 21. Cusped xi orbits in the instantaneous frozen rotating potential, superposed on the spatial gas distribution of 
model 110t2066 (left), and the corresponding traces in the £ — V diagram (right). The solid line orbit is obtained imposing 
reflection symmetry on the potential with respect to the vertical axis through the centre of mass, and the other orbits are based 
on the exact potential rotating about the centre of mass, the dashed and dotted orbits having single cusps at opposite ends. The 
crosses indicate the Lagrangian points Li, L2, L4 and L5, and the thin solid lines the corotation circle, based on the radius of 
the Lagrangian points Li and L2, and the longitudes of its tangent points. Note that in the particle representation used here, the 
£ — V plot looks different from the one in Fig. [l^ because the weights of the inverse squared distance relative to the observer are 
not taken into account. For example, the molecular ring is more difficult to distinguish without this correction. 



rather passes round of it and lands on the opposite dustlane 
closer to the nuclear ring, merging with it only at the next pas- 
sage or after repeating the whole cycle once more (Fig. ^ see 
also Fig. 10a in Fukuda et al. 1998). The farther out the gas 
on the dustlanes originates, the more it passes away from the 
nuclear ring. The upper and lower edges of the 180-pc parallel- 
ogram could represent such gas streams which precisely brush 
the nuclear ring/disc, loosing only little mass to it. The mass 
transfer would appear in the £ — V plot like vertical bridges be- 
tween the streams and the nuclear disc. An example of such 
bridges is detectable in the high resolution ^^CO data (e.g. 
Morris & Serabyn 1996) near i = 1.3° and for V between 100 
and 200 km s~^ . The brushing streams finally crash in the dust- 
lanes where they are abruptly decelerated. The right longitudi- 
nal edge of the 180-pc parallelogram, at £ —0.8°, could be a 
trace of this process. However it is not very clear why this edge 
is located at lower absolute longitude than the bright emission 
near the negative terminal velocity peak, corresponding to gas 
with about maximum velocity on the far-side dustlane, i.e. why 
the gas on this dustlane slows down before being striked by 
the positive velocity stream round the nuclear ring. Uchida et 
al. (1994b) have reported on a large scale shock front in the 
AFGL 5376 region, which is close to the southern end of the 
parallelogram upper edge, but which they consider as a support 
of the expanding molecular ring hypothesis. 



7. Geometrical constraints on the bar parameters 

Figure ^ shows how the ^ — F diagram of model 110' 12540 
changes when the viewing point of the observer is modified. 
Reducing the bar inclination angle ifo shrinks the structures 
longitudinally and amplifies the velocities. The observed knee 
of the 135-km5~^ arm near £ = — 5° is better reproduced, but 
the transition of the 3 -kpc arm to the connecting arm happens 
at too negative velocity and the connecting arm becomes too 
steep. Increasing the angle (^o lowers the forbidden velocities, 
moves the terminal velocity peaks further away from £ = and 
shifts the connecting arm closer to the northern tangent points 
of the two molecular ring branches. Increasing the distance of 
the observer obviously produces a longitudinal contraction, but 
without modifying the velocity of the structures near the centre. 
The diagram with the observer at = 8 moves the loop asso- 
ciated to the ^ < prolongation of the molecular ring structure 
out to the real tangent point of the Carina arm near I = —78°. 
These properties cannot be used to infer a robust inclination 
angle of the bar because they are based on one specific model 
and the gas flow is strongly time-dependent. 

Constraining the bar parameters by adjusting gas dynami- 
cal models to the observed CO and Hli — V diagrams is a very 
delicate task and may lead to unreliable results if the models 
are not sufficiently realistic. However, with our interpretation 
of the dominant features in these diagrams, it is possible to pro- 
vide geometrical constraints on the bar inclination angle and 
extension which do not depend on the details of the models. 
The principle of the method, illustrated in Fig. |3[ is to deter- 
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Fig. 22. Dependence of the ^ — F diagrams on the location of the observer relative to the bar in model 110't2540. The upper 
diagrams illustrate the effects of changing the bar inclination angle and the lower ones the effects of distance, represents the 
galactocentric distance of the observer in initial units, which is 9 in the adopted optimum model displayed in Fig. [iSl 



mine in the CO and HI data the longitudes ii and £2 where 
the 3-kpc and the 135-kms~^ arms intersect the axis shocks 
and to adjust in real space a major axis through the Galactic 
centre crossing the line of sights associated to these two direc- 
tions at galactocentric distances in the ratio q = 51/52. Simple 
trigonometrical considerations in the first and fourth Galactic 
quadrants respectively yield: 



sin \ii 



1, 2, 



Ro sin {cpo + ii) 
and isolating (fo and si between these two equations: 
1 + ^ 



tan(^o 
Ro 



cot £2 -\- Q cot ii ' 

sin^ ii -\- q sin^ £2 



(1 + ^)- 



sin^ {ii - £2) 



1/2 



(26-27) 



(28) 



(29) 



where q remains as a parametrisation of the asymmetry level 
between the lateral arms, which is not a priori known. In the 



ideal bisymmetric case q = I, but in reality q > 1 according to 
Sect. These formulae rest on the implicit assumption that 
the intersections of the lateral arms with the axis shocks and 
the Galactic centre are collinear. If this is wrong, then cpo will 
represent the angle between the line joining these intersections 
and the direction £ = 0. 

Extrapolating the connecting arm down to the 3-kpc arm in 
the £ — V observations leads to £1 ^ 14.5°, where a knee of 
the 3-kpc arm can be barely detected, whereas the 135-kms~^ 
arm meets the opposite axis shock at ^2 ^ —4.5° (Fig. ^. The 
resulting constraints are put together in Fig. Clearly, the 
distance si increases for smaller values of (fo, but no precise 
value can be given so far for the inclination angle. However, 
the CO data in Fig. [l| betray a second fainter far- side lateral 
arm which extends down to longitude £2 = —7° ±0.5° with 
lower forbidden velocities than the 135-kms~^ arm, as well 
as a quasi- vertical feature at about the same longitude proba- 
bly corresponding to gas from the same arm plunging towards 
the nuclear ring/disc after apocentre, on orbits parallel to the 
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Fig. 24. Graphical representation of the geometrical constraints 
on the bar parameters. The solid curves are derived from 
Eq. (26) for different values of , and the grey surfaces from 
Eq. (27) for different values of q, setting 52 = and assum- 
ing —4.75° < ^2 ^ —4.25° to show how uncertainties in £2 
propagate to q. 



Fig. 23. Basic geometry for the localisation of the intersections 
between the lateral arms (solid lines) and the axis shocks (long 
dashed lines). The nuclear ring (short dashed line) has been en- 
larged and the location of the Sun is indicated by the symbol. 



main axis shock. If this is correct, this arm must be much more 
symmetrical to the 3-kpc arm and the formulae ( ^8| ) and ( ^ ) 
can be applied to these two arms using q = 1. The results are 
LPo = 25° ± 4° and si = 3.15 kpc, and by the way q = 1.76, 
i.e. S2 = 1.8 kpc, for the 135-kms~^ arm. 

It remains to see how the arm intersections discussed here 
are related to the true bar parameters. Numerical simulations 
and analyses of observations in early-type barred galaxies in- 
dicate that the ratio between the bar semi-major axis and the 
corotation radius amounts to a/RL = 0.85 ± 0.15, and off- 
set dustlanes in such galaxies do not extend beyond the bar 
ends, i.e. Si < a (e.g. the review of Elmegreen 1996). Also, 
bisymmetric hydro simulations in rotating barred potentials 
with straight offset dustlanes generally have lateral arms in- 
tersecting the axis shocks very close to their outer ends and at 
most a few degrees ahead of the bar major axis. In the stan- 
dard model 001 of Athanassoula (1992), which has an inner 
Lindblad resonance but no looped xi orbits, s/Rl ~ 0.72. 
Adopting si/Rl = 0.8±0.1asa good compromise, our value 
of 5i would imply a corotation radius of 4.0 ± 0.5 kpc for the 
Milky Way. 

The situation in our non- symmetrised and time-dependent 
simulations, shown in Fig. is much more complicated how- 
ever. As mentioned in Sect. 5.2, the radius of the intersections 



between the lateral arms and the axis shocks increases with 
time. At the formation of a lateral arm, the associated intersec- 
tion leads the bar major axis, and as the arm moves outwards, it 



crosses this axis and becomes trailing. When both are exactly 
in phase, the intersection radius relative to corotation is 0.75 
on the average, compatible with the value adopted above, and 
compares very well to the apocentre radius of the cusped xi or- 
bit (see Fig. |2l|). Further out, the lateral arms rapidly dissolve in 
the spiral arms emanating from the very end of the axis shocks. 
The minimum value of the ratio Si/ Rl, which could depend 
on the size of the nuclear ring in the models, is about 0.4. Tak- 
ing this value as a lower limit for the 135-kms~^ arm leads to 
si/Rl = q ■ 82/ Rl ^ 0.7 for the 3-kpc arm, suggesting that 
this arm should meet the connecting arm close to the bar ma- 
jor axis and close to the apocentre of the cusped xi orbit. This 
is also confirmed by the fact that there is no obvious velocity 
gap in the observed i — V diagrams at the longitude where the 
3-kpc arm passes into the connecting arm. Our models suggest 
that the tangent points of the two molecular ring branches in the 
first Galactic quadrant are within corotation (see Fig. |2l|), con- 
trary to Englmaier & Gerhard's (1999) deduction, and therefore 
we rather advocate i^L > 4 kpc. 

Our value for the bar inclination angle agrees very well with 
many other determinations (e.g. Stanek et al. 1997; Nikolaev & 
Weinberg 1997), and in particular with the recent result (po = 
20° - 25° of Englmaier & Gerhard (1999), who derived the gas 
flow in the potential of the deprojected COBE/DIRBE L-band 
luminosity distribution (Binney et al. 1997). The bisymmet- 
ric hydro simulations done by Weiner & Sellwood (1999) and 
matched to the observed HI terminal velocity curves favour 
slightly larger values of this angle, i.e. (fo = 35° ± 5°, and 
support a corotation radius between 4 and 6 kpc (for Ro = 
8.5 kpc). These authors consider an angle below 25° unlikely, 
but it should be noted that gas flow fluctuations like those in 
our self-consistent simulations can increase the£ — V area cov- 
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Fig. 25. Location of the intersections between the lateral arms and the axis shocks relative to the bar in simulation 110. Left top: 
distance of these intersections relative to the offcentred density centre and normalised by the corotation radius (based on the 
Lagrangian points Li and L2). Left bottom: phase of the intersections relative to the major axis of the bar, defined positive when 
the formers lead the latter; the dotted line indicates the locus where both are exactly in phase. Right: correlation between the 
ordinates of the two previous plots. The solid and dashed lines refer to arm intersections on the same side of the bar, and the bold 
lines to intersections involving the most contrasted lateral arms. 



ered by forbidden velocities and hence substantially bias their 
method. 

It is also possible from our data interpretation to give a 
crude estimate of the bar pattern speed. If the 3-kpc arm indeed 
encounters the connecting arm near the apocentre of the cusped 
xi orbit, then the gas at the intersection of the two arms should 
corotate with the bar figure. Evaluating the radial velocity AV 
of this gas relative to the LSR yields its rotation velocity with 
respect to the Galactic centre: 

AV + Ksin^i 



(30) 



and resorting to Eq. (26): 

s\ itoSmi;i 

where Vo is the circular velocity of the LSR and Qo = Vo/Ro. 
The bar inclination angle and si have vanished in this last for- 
mula and thus the method is independent of them. The value 
of /SV is hard to determine in the CO and Y{\ I — V plots 
because of overlayed emission from the disc and in particu- 
lar from the molecular ring, but a reasonable range is = 
15 - 35 kms-^ With i?o = 8 kpc, K = 220 kms"^ and 
ii = 14.5° as before, one gets l]p = 35 — 45 km s~^ kpc~^, in 



fair agreement with the 50 km s~^ kpc~^ of our gas simulations 
(Fig. [1^) when rescaled to = 4.25 kpc. Changing the value 
of £1 or by 1° and 10 kms~^ modifies the result only by 
~ 2% and 3% respectively. However, AV may be affected by a 
radial motion of the LSR and/or oscillations of the bar density 
centre. Longitude-velocity maps of very dense gas (like CS) 
could help to localise more precisely the transition from the 
3-kpc arm to the connecting arm. Note that the method does 
not apply to the very asymmetric model 110' 12540, as can be 
checked from Fig. ^ because the gas at the transition is not at 
rest in the frame of the bar. 

8. Conclusion 

This paper presents the first fully self-consistent gas flow mod- 
els of the inner Milky Way, obtained by incorporating an SPH 
component with 150 000 particles in realistic symmetry-free 
barred 3D A^-body simulations with nearly 4x10^ total number 
of particles. The axisymmetric initial conditions of the simula- 
tions are chosen to evolve spontaneously in a barred configu- 
ration compatible with the COBE/DIRBE K-band constraints, 
according to criteria based on a set of lower resolution pure 
stellar dynamical simulations realised in a precedent paper. The 
stellar bar rotates with its intrinsic natural pattern speed and the 
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gas component, gently released from its axisymmetric config- 
uration after the formation of the bar, freely interacts with the 
live stellar arms. 

The density centre of the stellar bar becomes unstable and 
oscillates around the centre of mass with an amplitude of sev- 
eral 100 pc and a frequency of 20 — 30 km s~^ kpc~^. The gas 
flow, contrary to other bisymmetric hydro simulations in fixed 
rotating barred potentials, never reaches a quasi- steady state 
and can delineate strong asymmetries in its spiral structure, not 
necessarily induced by the lopsided stellar distribution. 

Some models selected from the simulations account for 
many features seen in the HI and CO longitude- velocity dis- 
tributions within the Galactic bar and surrounding regions, and 
provide a very powerful guide to understand the inner structure 
of the Milky Way. The Galactic bar is inclined by 25° ± 4° rel- 
ative to the £ = line, has a corotation radius of 4.0 — 4.5 kpc, 
a related pattern speed ~ 50 kms~^ kpc~^, and a face-on 
axis ratio 6/a 0.6. As in most early-type barred spirals, offset 
dustlanes are leading the bar major axis. Their gaseous traces 
in the observed i — V diagrams correspond to the connecting 
arm (near- side branch) and another feature near ^ = — 4° (far- 
side branch). These dustlanes are the loci of strong shearing 
shocks with velocity changes up to 200 km s~^ across, and are 
located closer to the bar major axis than the leading edges of 
the cusped xi orbit. The peaks in the terminal velocity curves 
at ^ ^ ±2.5° are produced by the post- shock gas and not by 
the £ — V trace of the latter orbit, which has velocity peaks 
at larger absolute longitude and with lower velocity amplitude. 
The near-side branch of the dustlanes lies below the Galactic 
plane (b < 0) and the other branch, which is seen almost end- 
on, above it. Their maximum departure fromthe plane amounts 
to more than 100 pc. 

The 3 -kpc and 135-kms~^ arms are the inner prolongations 
of spiral arms in the disc, each passing close to one bar end 
and joining by a large bow around the nuclear ring/disc the 
dustlane in the other side of the bar, at galactocentric distances 
R ^ 3.2 kpc and 1.8 kpc respectively. The 135-kms~^ arm is 
associated with larger forbidden velocities because of its deeper 
location in the central potential well. 

Velocity-elongated features in the observed CO ^ — F di- 
agrams "below" the connecting arm, i.e. at similar longitudes 
but lower velocities than this arm, are interpreted as gas lumps 
which are just about to cross the shock front layer of the near- 
side dustlane. A robust example of such a feature is given by 
the vertical feature near £ = 5.5°. Bania's clump 2 complex is a 
more controversial candidate because its mean Galactic latitude 
differs by roughly half a degree from that of the connecting arm 
at same longitude and because of its substantial mass. Finally, 
Bania's clump 1 complex is the part of the 135-kms~^ arm 
which strikes the far- side dustlane with a considerable incident 
velocity of order 100 kms~^. A fraction of its gas is directly 
absorbed by the dustlane, producing the apparent positive ve- 
locity trace of the dustlane, another fraction is gliding outwards 
along the dustlane, but the bulk of the clump is probably not in- 
tercepted by the dustlane. 



The interpretation of the main observed ^— ^ features given 
in this paper are by far the most plausible that can be deduced 
from our simulations. However, these simulations may still sig- 
nificantly depart from reality and hence the unicity of our inter- 
pretation cannot be guaranteed. For instance, if a large amount 
of gas is concentrated in the centre, a nuclear fast rotating bar 
may temporarily form and lead to a different understanding of 
the very complex innermost gas dynamics. 

Some mpeg movies of the gas flow in simulation 110, 
including live £ — V di agrams, are available on the web at 



|http://obswww.unige.chy ~fux. 
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